for u = [6,58]
Import Specific Image
I=imread(char(strcat({'0101_baseline_anterior copy'}, {' '}, {num2str(u)}, {'.jpg'})));
imshow(I)
Crop & Convert to Blackwhite
J=I;
%J=imcrop(I,[0.5 40.5 613 397]);
J=rgb2gray(J);
imshow(J)
Select the degree celcius symbol, then detect where it exists.
I2=imread('0101_baseline_anterior copy 6.jpg');
J2=rgb2gray(I2);
J2(8:32,149:174)=J2(8:32,149:174)*1.16;
J_test=J2>253;
J_test(35:450,1:590)=0;
cc = bwconncomp(J_test);
L = labelmatrix(cc);
bw2 = bwdist(J_test) <= 25;
L4 = labelmatrix(bwconncomp(bw2));
L4(~J_test) = 0;
%L4(:,[1:152 171:640]) = 0;
%L4(40:480,152:171) = 0;
%L4=logical(L4);
L4=logical(L4);
imshow(L4)
Loop to reverse-shade shaded regions (should also test for if regions exist and are equivelent to the ones defined here...)
if sum(sum(immultiply(J,L4)))/sum(sum(L4))>200
x1=[6, 149, 6, 556, 556];
x2=[301, 174, 141, 635, 635];
y1=[450, 6, 6, 441, 6];
y2=[474, 34, 38, 475, 40];
else
x1=[6, 13, 556, 556];
x2=[301, 38, 635, 635];
y1=[450, 6, 441, 6];
y2=[474, 34, 475, 40];
end
for k=1:length(x1)
if length(x1)==5
if k==1 || k==2 || k==3
s=2;
else s=1;
end
else
if k==1 || k==2
s=2;
else s=1;
end
end
for i = (x1(k)+2):(x2(k)-2)
psi{k}(i,1)=mean(J((y1(k)-1-s):(y1(k)-1),i));
psi{k}(i,2)=mean(J((y2(k)+1):(y2(k)+1+s),i));
psi{k}(i,3)=mean(J((y1(k)):(y1(k)+s),i));
psi{k}(i,4)=mean(J((y2(k)-s):(y2(k)),i));
for j = y1(k):y2(k)
J(j,i)=J(j,i)*((psi{k}(i,1)/psi{k}(i,3))*((y2(k)-j)/(y2(k)-y1(k)))+...
(psi{k}(i,2)/psi{k}(i,4))*((j-y1(k))/(y2(k)-y1(k))));
end
end
for i = [(x1(k)+1),(x2(k)-1)]
psi{k}(i,1)=mean(J((y1(k)-1-s):(y1(k)-1),i));
psi{k}(i,2)=mean(J((y2(k)+1):(y2(k)+1+s),i));
psi{k}(i,3)=mean(J((y1(k)+1):(y1(k)+1+s),i));
psi{k}(i,4)=mean(J((y2(k)-1-s):(y2(k)-1),i));
for j = (y1(k)+1):(y2(k)-1)
J(j,i)=J(j,i)*((psi{k}(i,1)/psi{k}(i,3))*((y2(k)-j)/(y2(k)-y1(k)))+...
(psi{k}(i,2)/psi{k}(i,4))*((j-y1(k))/(y2(k)-y1(k))));
end
end
for i = [(x1(k)),(x2(k))]
psi{k}(i,1)=mean(J((y1(k)-1-s):(y1(k)-1),i));
psi{k}(i,2)=mean(J((y2(k)+1):(y2(k)+1+s),i));
psi{k}(i,3)=mean(J((y1(k)+2):(y1(k)+2+s),i));
psi{k}(i,4)=mean(J((y2(k)-2-s):(y2(k)-2),i));
for j = (y1(k)+2):(y2(k)-2)
J(j,i)=J(j,i)*((psi{k}(i,1)/psi{k}(i,3))*((y2(k)-j)/(y2(k)-y1(k)))+...
(psi{k}(i,2)/psi{k}(i,4))*((j-y1(k))/(y2(k)-y1(k))));
end
end
end
imshow(J)
Find area of text & connect nearby regions.
J_test=J>253;
J_test(35:450,1:590)=0;
cc = bwconncomp(J_test);
L = labelmatrix(cc);
rgb = label2rgb(L, 'jet', [.7 .7 .7], 'shuffle');
imshow(rgb)
bw2 = bwdist(J_test) <= 25;
imshow(bw2)
L2 = labelmatrix(bwconncomp(bw2));
rgb2 = label2rgb(L2, 'jet', [.7 .7 .7], 'shuffle');
imshow(rgb2)
L2(~J_test) = 0;
imshow(label2rgb(L2, 'jet', [.7 .7 .7], 'shuffle'))
Remove text by dilating (outward) the just text, then inward interpolation
L3=imbinarize(L2*255);
se = strel('rectangle',[4 4]);
dilatedL3 = imdilate(L3,se);
dilatedL3 = imbinarize(dilatedL3*255);
imshow(dilatedL3)
JJ = immultiply(J,~dilatedL3);
JJ = double(JJ);
JJ(JJ==0)=NaN;
JJ=inpaint_nans(JJ,5);
JJ=uint8(JJ);
imshow(JJ)
Remove verticle bar on the right of the image
L4=zeros(480,640);
L4(42:439,615:636)=1;
JJ = immultiply(JJ,~L4);
JJ = double(JJ);
JJ(JJ==0)=NaN;
JJ=inpaint_nans(JJ);
JJ=uint8(JJ);
imshow(JJ)
centreregion = zeros(480,640);
centreregion(172:280,167:360)=1;
cross = JJ<50;
cross = immultiply(cross, centreregion);
if sum(sum(cross)) >= 10
se = strel('rectangle',[4 4]);
dilatedcross = imdilate(cross,se);
JJ = immultiply(JJ,~dilatedcross);
JJ = double(JJ);
JJ(JJ==0)=NaN;
JJ=inpaint_nans(JJ, 4);
JJ=uint8(JJ);
end
imshow(JJ)
Test if are present in J, then remove image's background:
First, compute x- & y-Derivatives
[J_Gx,J_Gy] = imgradientxy(JJ);
J_Gx=abs(J_Gx);
imshow(J_Gx)
Perform first Otsu's Test (does ?) to remove background & disply output beside original image
level = graythresh(JJ);
BW = imbinarize(JJ,level);
BW=imfill(BW,'holes');
nBW=imfill(~BW,'holes');
BW=immultiply(BW,~nBW);
imshowpair(JJ,BW,'montage')
Test if , and , then if true, perform corrective filling-in procedure to properly segment the background:
krn=[1 -1];
dimBW=size(BW);
if ~(mean(mean(J_Gx(1:5,:)))>mean(mean(J_Gx)))...
&& ((sum(sum(~imcrop(BW,[238.5 335.5 114 52])))>0 ...
|| sum(sum(~imcrop(BW,[152.5 121.5 72 72])))>0 ...
||sum(sum(~imcrop(BW,[402.5 141.5 72 72])))>0))
for i = 1:dimBW(1)
changes_first{i}=conv(krn, BW(i,:));
changes1_first{i}=find(changes_first{i}==1 | changes_first{i}==-1);
end
BW(1,changes1_first{1}(1):changes1_first{1}(length(changes1_first{1})-1))=1;
BW(dimBW(1),changes1_first{dimBW(1)}(1):...
changes1_first{dimBW(1)}(length(changes1_first{dimBW(1)})-1))=1;
BW = imfill(BW, 'holes');
imshowpair(JJ,BW,'montage')
end
If , and was true, test if previous code corrected the problems:
if ~((sum(sum(~imcrop(BW,[238.5 335.5 114 52])))==0 ...
&& sum(sum(~imcrop(BW,[152.5 121.5 72 72])))==0 ...
&& sum(sum(~imcrop(BW,[402.5 141.5 72 72])))==0))
fprintf('WARNING: Otsu Test Failure - Complete Dependence on Alternative to Otsu');
end
Automatically detect if Otsu's test failed to remove the proper background, and if so apply 2nd method for removing the background [so far dependent on mu & sigma in:
edge(JJ, 'Canny', mu ,sigma) , and dependent on 4 distinct background regions - need to generalize]
if (sum(sum(~imcrop(BW,[238.5 335.5 114 52])))>0 ...
|| sum(sum(~imcrop(BW,[152.5 121.5 72 72])))>0 ...
||sum(sum(~imcrop(BW,[402.5 141.5 72 72])))>0)
BW1 = edge(JJ,'Canny',.7,3);
s = regionprops(BW1,'centroid');
centroids = cat(1, s.Centroid);
imshow(BW1)
hold on
plot(centroids(:,1),centroids(:,2), 'g*')
hold off
s1 = regionprops(BW1,'Extrema');
s3 = regionprops(BW1,'BoundingBox');
s1 = struct2cell(s1);
s3 = struct2cell(s3);
dim=size(BW1);
BW2=BW1;
dim3=size(s3);
if dim3(2)==3 && centroids(2,2)<mean(centroids(:,2))
%gamma1
BW2(floor(s3{1}(2)):dim(1),1)=1;
BW2(dim(1),1:(s3{1}(3)))=1;
%gamma4
BW2(floor(s3{3}(2)):dim(1),dim(2))=1;
BW2(dim(1),floor(s3{3}(1)):dim(2))=1;
%gamma2
BW2(1,floor(s3{2}(1)):(floor(s3{2}(1))+s3{2}(3))+1)=1;
end
if dim3(2)==4 && centroids(2,2)<mean(centroids(:,2))...
&& centroids(3,2)<mean(centroids(:,2))
%gamma1
BW2(floor(s3{1}(2)):dim(1),1)=1;
BW2(dim(1),1:(s3{1}(3)))=1;
%gamma4
BW2(floor(s3{4}(2)):dim(1),dim(2))=1;
BW2(dim(1),floor(s3{4}(1)):dim(2))=1;
%gamma2
BW2(1,floor(s3{2}(1)):(floor(s3{2}(1))+s3{2}(3))+1)=1;
%gamma3
BW2(1,floor(s3{3}(1)):(floor(s3{3}(1))+s3{3}(3))+1)=1;
end
if (dim3(2)==3 && centroids(2,2)<mean(centroids(:,2))...
&& centroids(3,2)<mean(centroids(:,2)))||...
(dim3(2)==4 && centroids(2,2)<mean(centroids(:,2)))
BW2 = imfill(BW2, 'holes');
BW2 = ~BW2;
imshow(BW2)
BW=BW2;
end
if ~(dim3(2)==3 || dim3(2)==4)
fprintf('FATAL WARNING: x=Regions(BW1) not in {3,4}');
end
if dim3(2)==3 && centroids(2,2)>mean(centroids(:,2))
fprintf('FATAL WARNING: x_2=Region_2(BW1) not < mean(Regions(BW1))');
end
if dim3(2)==4 && (centroids(2,2)>mean(centroids(:,2))||...
centroids(3,2)>mean(centroids(:,2)))
fprintf('FATAL WARNING: x_2=Region_2(BW1) or x_3=Region_3(BW1) not < mean(Regions(BW1))');
end
end
Restrict Calculations to the proper subset of the Image:
JJJ=immultiply(JJ,BW);
JJJ(JJJ==0)=NaN;
imshow(JJJ)
dimIJK=size(JJJ);
IJK=logical(JJJ);
imshow(IJK)
Reduce Cropping (x-direction dependent) by 3 pixels
krn=[1 -1];
for i = 1:dimIJK(1)
changes{i}=conv(krn, IJK(i,:));
changes1{i}=find(changes{i}==1 | changes{i}==-1);
end
for i = 1:dimIJK(1)
IJK(i, max(changes1{i}-1,1))=0;
IJK(i, max(changes1{i}-2,1))=0;
IJK(i, max(changes1{i}-3,1))=0;
IJK(i, max(changes1{i}-4,1))=0;
IJK(i, min(changes1{i}+1,dimIJK(2)))=0;
IJK(i, min(changes1{i}+2,dimIJK(2)))=0;
IJK(i, min(changes1{i}+3,dimIJK(2)))=0;
IJK(i, min(changes1{i}+4,dimIJK(2)))=0;
IJK(i, min(changes1{i}, dimIJK(2)))=0;
end
JJ=immultiply(JJ,uint8(IJK));
JJ(JJ==0 | isnan(JJ))=NaN;
imshow(JJ)
Compute x- & y-Derivatives
[Gx,Gy] = imgradientxy(JJ);
Gxy=Gx+Gy;
Gy=abs(Gy);
Gx=abs(Gx);
Gxy=abs(Gxy);
Gxy=imgaussfilt(imgaussfilt(Gxy));
Gy=imgaussfilt(imgaussfilt(Gy));
Gx=imgaussfilt(imgaussfilt(Gx));
Gx2=Gx;
Gx2(:,[1:224,416:640])=0;
imshow(Gx)
Find the most most significant y-derivative areas
First Guess: .3*nanstd(nanstd(double(Gy)))^.5*range(range(Gy))
totest=Gy+.1*Gx+Gx2+.6*Gxy;
IJ=totest>=max(max(Gy))-.92*floor(range(range(Gy)));
dimIJ=size(IJ);
IJ(1:floor(0.48*dimIJ(1)),1:dimIJ(2))=0;
imshow(IJ)
Create the Hough transform using the binary image.
[H,T,R] = hough(IJ);
clear step
for i = 1:180
step{i} = find(logical(H(:,i)));
H(step{i}(1):(step{i}(1)+floor(.2*(step{i}(length(step{i}))-step{i}(1)))),i)=0;
H((step{i}(1)+floor(.8*(step{i}(length(step{i}))-step{i}(1)))):step{i}(length(step{i})),i)=0;
end
imshow(H(:,1:180),[],'XData',T,'YData',R,...
'InitialMagnification','fit');
xlabel('\theta'), ylabel('\rho');
axis on, axis normal, hold on;
Find peaks in the Hough transform of the image.
clear P1 P2 P
for q = 1:2
P1{1} = houghpeaks(H(:,1:75),4);
P2{1} = houghpeaks(H(:,106:180),4);
P1{2} = houghpeaks(H(:,1:6),3);
P2{2} = houghpeaks(H(:,175:180),3);
P2{1}(:,2) = 90+P2{1}(:,2);
P2{2}(:,2) = 174+P2{2}(:,2);
P = vertcat(P1{q},P2{q});
x = T(P(:,2)); y = R(P(:,1));
plot(x,y,'s','color','white');
hold off
Find lines and plot them.
clear lines euclid euclid2 euclid3
if q==1
lines = houghlines(IJ,T,R,P,'FillGap',6, 'MinLength', 8);
else
lines = houghlines(IJ,T,R,P,'FillGap',3, 'MinLength', 8);
end
euclid1 = 1:length(lines);
for i = 1:length(lines)
if lines(i).theta<0
euclid1(i)=1/2*lines(i).point1(1)+1/2*lines(i).point1(1)>305;
else
euclid1(i)= 1/2*lines(i).point1(1)+1/2*lines(i).point1(1)<335;
end
end
euclid = find(euclid1);
for k = 1:length(lines)
if euclid1(k)==1
euclid(k) = pdist([lines(k).point1; lines(k).point2],'euclidean');
else
euclid(k) = 0;
end
end
euclidn = euclid.*([lines.theta]<0);
euclidp = euclid.*([lines.theta]>=0);
euclid2n = sort(euclidn);
euclid2p = sort(euclidp);
if q == 1
pts = 6; % number of lines you want, must also be even
else
pts = 2; % as a standard, pts | q=2 will be 2
end
euclid3n = zeros(1,length(lines));
for k = find([lines.theta]<0)
if euclid(k) > euclid2n(length(lines)-pts/2)
euclid3n(k)=1;
else
euclid3n(k)=0;
end
end
euclid3p = zeros(1,length(lines));
for k = find([lines.theta]>=0)
if euclid(k) > euclid2p(length(lines)-pts/2)
euclid3p(k)=1;
else
euclid3n(k)=0;
end
end
euclidall=logical(euclid3n+euclid3p);
lines = lines(euclidall);
figure, imshow(IJ), hold on
max_len = 0;
for k = 1:length(lines)
xy = [lines(k).point1; lines(k).point2];
plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','green');
% Plot beginnings and ends of lines
plot(xy(1,1),xy(1,2),'x','LineWidth',2,'Color','yellow');
plot(xy(2,1),xy(2,2),'x','LineWidth',2,'Color','red');
% Determine the endpoints of the longest line segment
len = norm(lines(k).point1 - lines(k).point2);
if ( len > max_len)
max_len = len;
xy_long = xy;
end
end
hold off
clear s
for k = 1:size(lines, 2)
x1 = lines(k).point1(1);
x2 = lines(k).point2(1);
y1 = lines(k).point1(2);
y2 = lines(k).point2(2);
if x1 < x2
s{k}(:,1)=x1:x2;
for i = 1:(x2-x1+1)
s{k}(i,2)=floor((y2-y1)/(x2-x1)*(s{k}(i,1)-x1)+y1);
s{k}(i,3)=ceil((y2-y1)/(x2-x1)*(s{k}(i,1)-x1)+y1);
end
else
s{k}(:,1)=x1:-1:x2;
for i = 1:(x1-x2+1)
s{k}(i,2)=floor((y2-y1)/(x2-x1)*(s{k}(i,1)-x1)+y1);
s{k}(i,3)=ceil((y2-y1)/(x2-x1)*(s{k}(i,1)-x1)+y1);
end
end
end
%x1 = min(lines(k).point1(1),lines(k).point2(1));
%x2 = max(lines(k).point1(1),lines(k).point2(1));
%y1 = min(lines(k).point1(2),lines(k).point2(2));
%y2 = max(lines(k).point1(2),lines(k).point2(2));
SI{q} = zeros(480, 640);
for k = 1:size(lines, 2)
SI{q}(480*(s{k}(:,1)-1)+s{k}(:,2))=1;
SI{q}(480*(s{k}(:,1)-1)+s{k}(:,3))=1;
end
SI{q} = logical(SI{q});
imshow(SI{q})
se = strel('sphere',4);
dilatedSI{q} = imdilate(SI{q},se);
imshow(dilatedSI{q})
end
dSI = or(dilatedSI{1}, dilatedSI{2});
IJ2 = immultiply(IJ,dSI);
IJ21=bwareafilt(IJ2,2);
imshow(IJ2)
s2 = regionprops(IJ21,'centroid');
centroids = cat(1, s2.Centroid);
imshow(IJ2)
middle = ([sum(centroids(:,1))/2,sum(centroids(:,2))/2]);
hold on
plot(centroids(:,1),centroids(:,2), 'rs', 'MarkerSize', 15)
plot(middle(1), middle(2),'gx', 'MarkerSize', 15)
hold off
leftregion = imcrop(IJ2,[0 0 middle(1) 480]);
rightregion = imcrop(IJ2,[middle(1) 0 640 480]);
Fit polynomial line to the left area
ii=find(leftregion);
[yy,xx]=ind2sub(size(IJ2),ii);
p1 = polyfit(xx,yy,2);
Right area
ii=find(rightregion);
[yy,xx]=ind2sub(size(IJ2),ii);
p2 = polyfit(xx,yy,2);
Plots of polynomial fits & saves it to cd
imshow(JJ)
x1 = 0:640;
y1 = polyval(p1,x1);
x2 = (-middle(1)):(640-middle(1));
y2 = polyval(p2,x2);
hold on
plot(x1,y1)
plot(x2+middle(1),y2)
saveas(gcf,strcat('poly',num2str(u)), 'jpg')
hold off
end